The goal of this notebook is to phase more clusters of variants using the same approach we used to define the ‘genotypes’ Genome-1 and Genome-2 in the previous notebook. The output of this notebook is the list of variants annotated by whether they are part of any haplotypes/clusters.

The inputs for this analysis are (1) the filtered and combined variants annotated by major genotypes from 02-determine-main-genotypes.Rmd and (2) annotations of the positions of the main MeV open reading frames.

## ==== File paths input data ==== ##

if (exists("snakemake")) {
  
  # Variant data annotated with major genotypes
  labeled.data = snakemake@params[['incsv']]
  
  # Annotations for the standard Measles ORFs  
  annotations.filepath = snakemake@params[['annotations']]
  
  # Path to the ancestral mutations
  ancestral.data = snakemake@params[['ancestral']]

} else {
  
  # Variant data annotated with major genotypes
  labeled.data = "../../results/variants/genotyped_variants.csv"
  
  # Annotations for the standard Measles ORFs  
  annotations.filepath = "../../config/annotations.csv"
  
  # Path to the ancestral mutations
  ancestral.data = "../../config/ref/annotated_SSPE_consensus_snps.csv"

}
## ==== File paths output data ==== ##

# Path to save the results
if (exists("snakemake")) {
  
  # Variants annotated by new haplotype clusters
  output.path = snakemake@params[['outcsv']]
  
  # Path to save figures
  figure.dir = paste0(snakemake@params[['figures']], "/")

} else {
  
  # Variants annotated by new haplotype clusters
  output.path = "../../results/variants/clustered_variants.csv"
  
  # Path to save figures
  figure.dir = "../../results/figures/"
  
}

# Make the figure directory if it doesn't exist
if (!file.exists(figure.dir)) {
  dir.create(figure.dir, recursive = TRUE)
  print("Directory created.")
} else {
  print("Directory already exists.")
}
## [1] "Directory already exists."

Identify Subclonal Haplotypes

Since the goal of this notebook is use the same technique that I used to identify Genome-1, Genome-2, and Genome-1-1, but for mutations that aren’t necessarily found in every single tissue, I need to account for missing mutations. To do this, I’ll spread the data frame of frequencies and assign an allele frequency of 0 if a variant is missing from a tissue. This will allow us to cluster all of the mutations, even if they weren’t observed in a tissue sample.

I’ll write the approach I used in 02-determine-main-genotypes.Rmd into a reusable function. I can provide a list of SNPs, the full data frame, and the number of clusters I’m targeting.

cluster.snps <- function(list.of.snps, snp.df, n.clusters) {
  
  # SNP as column and Allele frequency as row
  frequency.by.snp = snp.df %>% 
    filter(SNP %in% list.of.snps) %>% 
    select(AF, SNP, Tissue) %>% 
    pivot_wider(names_from = SNP, values_from = AF, values_fill = 0) %>% 
    select(!Tissue) 
  
  # Calculate R between every pair of columns while handling NAs 
  snp.correlation = cor(frequency.by.snp, use="pairwise.complete.obs")
  
  # Convert to distance (positive corr is close to 0)
  snp.dist = as.dist(1 - snp.correlation)
  
  # Create k-medoids clustering with n clusters
  snp.kmedoids = pam(snp.dist, n.clusters) 
  kclusters = snp.kmedoids$cluster
  
  # Convert to a data frame
  kclusters.df = data.frame(SNP = names(kclusters), cluster = kclusters)

  # Assign the clusters to the original data frame
  kmedoids.SNPs = snp.df %>% 
    filter(SNP %in% list.of.snps) %>% 
    left_join(., kclusters.df, by = "SNP") %>% 
    mutate(cluster = if_else(is.na(cluster), "no cluster", as.character(cluster)))
  
  # Sort the cluster names by mean frequency 
  cluster.order = kmedoids.SNPs %>% 
    group_by(cluster) %>% 
    summarize(AF = mean(AF)) %>% 
    arrange(-AF) %>% 
    mutate(order = 1:n.clusters) %>% 
    select(cluster, order)
  kmedoids.SNPs = kmedoids.SNPs %>% 
    left_join(., cluster.order, by = "cluster") %>% 
    select(!cluster) %>% 
    rename(cluster = order) %>% 
    mutate(cluster = as.character(cluster))
  
  # Add the number of SNPs per cluster
  n.clusters.per.snp = kmedoids.SNPs %>% 
    select(SNP, cluster) %>%
    distinct() %>% 
    group_by(cluster) %>% 
    count() %>% 
    mutate(cluster_size = paste0(cluster, " (", n, " snps in cluster)"))
  
  return(left_join(kmedoids.SNPs, n.clusters.per.snp, by = "cluster"))
  
}

Although this approach works reasonably well, it’s very susceptible to noise. This means that you can’t just provide all the variants at once and expect the algorithm to partition them into correct clusters. Instead, the best approach I found was to break variants down into bins based on their frequency and cluster these separately. The higher average frequency, the less noise, and the better the algorithm is at partitioning the variants into meaningful clusters.

>= 25% Allele Frequency Mutations

First, I’ll start with a ‘high-frequency’ bin. This includes all variants that reach greater than or equal to 25% in at least one tissue.

Some clusters have only a single mutation in them. Clearly, some of these belong to Genome-1, Genome-2, and Genome-1-1, but are missing from a single tissue, perhaps due to low coverage in that particular tissue.

Cluster 5 (L-A15084G) is a mutation in Genome-1-1 that was too low frequency to be called in the Frontal Cortex 2 sample. Cluster 6 (M-C4502T) is mutation that’s part of Genome-2 but was too low coverage in the Cerebellum. Finally, Cluster 9 (M-C4573T) is a mutation in Genome-1 that was also too low coverage in the Cerebellum. Generally, the Cerebellum has the lowest coverage of any tissue, so this isn’t surprising. All-in-all, this adds three mutations to already existing haplotypes. I already addressed F-C7036T, N-G810A, and H-T7293C in the identify-backgrounds.Rmd notebook. Finally, M-C4532T might also belong as part of the reference?

I’ll annotate these variants in the main data frame.

# Clusters that are clearly genome-1 
genome.1.missing.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(4)) %>% 
  pull(SNP) %>% 
  unique(.)

# Clusters that are clearly genome-1-1
genome.1.1.missing.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(5)) %>% 
  pull(SNP) %>% 
  unique(.)

# Clusters that are clearly genome-2 
genome.2.missing.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(6)) %>% 
  pull(SNP) %>% 
  unique(.)

# Clusters that break the rules 
maybe.in.both.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(3, 7, 2, 1)) %>% 
  pull(SNP) %>% 
  unique(.)


# Annotate the labeled data frame with what's possible here
updated.haplotype.df = labeled.df %>%
  mutate(Background = case_when(SNP %in% genome.1.missing.snps ~ "genome-1",
                                SNP %in% genome.1.1.missing.snps ~ "genome-1",
                                SNP %in% genome.2.missing.snps ~ "genome-2",
                                SNP %in% maybe.in.both.snps ~ "both",
                                TRUE ~ Background)) %>% 
  mutate(Haplotype = case_when(SNP %in% genome.1.missing.snps ~ "genome-1",
                               SNP %in% genome.1.1.missing.snps ~ "genome-1-1",
                               SNP %in% genome.2.missing.snps ~ "genome-2",
                               SNP %in% maybe.in.both.snps ~ "both",
                               TRUE ~ Haplotype))

Now, lets look at the remainder of the clusters of SNPs. These are SNPs that actually clustered with other variants.

All of these look like reasonable clusters, so I’ll annotate the variant data frame with these new clusters (subclonal haplotypes). I’ll also try to guess if these are on the background of Genome-1 or Genome-2 if it’s possible.

cluster.1.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(8)) %>% 
  pull(SNP) %>% 
  unique(.)

# The missing SNPs in Cerebellum are due to coverage
cluster.2.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(11, 14)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.3.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(10)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.4.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(13)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.5.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(15)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.6.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(16)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.7.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(17)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.8.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(18)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.9.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(19)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.10.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(20)) %>% 
  pull(SNP) %>% 
  unique(.)
  
# Annotate the labeled data frame with what's possible here
updated.haplotype.df = updated.haplotype.df %>%
  mutate(Haplotype = case_when(SNP %in% cluster.1.snps ~ "cluster 1",
                               SNP %in% cluster.2.snps ~"cluster 2",
                               SNP %in% cluster.3.snps ~"cluster 3",
                               SNP %in% cluster.4.snps ~ "cluster 4",
                               SNP %in% cluster.5.snps ~ "cluster 5",
                               SNP %in% cluster.6.snps ~ "cluster 6",
                               SNP %in% cluster.7.snps ~ "cluster 7",
                               SNP %in% cluster.8.snps ~ "cluster 8",
                               SNP %in% cluster.9.snps ~ "cluster 9",
                               SNP %in% cluster.10.snps ~ "cluster 10",
                               TRUE ~ Haplotype))

Haplotypes >= 5% and < 25%

Now, I’ll move onto the next bin. These are mutations that reach at least 5% in one tissue but weren’t part of the previous analysis. Most of these will be challenging to haplotype due to noise. I’ll only try to specify the clusters that seem the most clear cut.

Most of these aren’t very obviously clusters that represent true haplotypes. They’re more likely just an artifact of noise. I took the 8 that looked the most promising. It’s probably possible to do this more systematically using the internal variation in a cluster, but the results will be the same and there aren’t enough clusters to justify the need for this approach.

cluster.11.snps = twentyfive.percent.or.fewer.clusters.df %>% 
  filter(cluster %in% c(15)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.12.snps = twentyfive.percent.or.fewer.clusters.df %>% 
  filter(cluster %in% c(4)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.13.snps = twentyfive.percent.or.fewer.clusters.df %>% 
  filter(cluster %in% c(14)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.14.snps = twentyfive.percent.or.fewer.clusters.df %>% 
  filter(cluster %in% c(12)) %>% 
  pull(SNP) %>% 
  unique(.)

# Annotate the labeled data frame with what's possible here
updated.haplotype.df = updated.haplotype.df %>%
  mutate(Haplotype = case_when(SNP %in% cluster.11.snps ~ "cluster 11",
                               SNP %in% cluster.12.snps ~ "cluster 12",
                               SNP %in% cluster.13.snps ~ "cluster 13",
                               SNP %in% cluster.14.snps ~ "cluster 14",
                               TRUE ~ Haplotype))

What’s left to cluster?

We were able to cluster a reasonable number of the high-frequency mutations (>= 5% frequency). How many mutations haven’t been haplotyped yet?

## [1] "There are about 471 that are still totally un-haplotyped."
## [1] "There are 60 that have been clustered."
## [1] "There are still about 262 that can reasonably be haplotyped (Above 5% in more at least one sample)."

All Current Haplotype Clusters

Here are the frequency of all current haplotypes in the 15 tissue samples.

Looking closely at the haplotypes, there is one amendment to be made. Cluster 10 and 14 are likely the same cluster. I think these weren’t clustered together because the mean frequency of these mutations is split between the frequency bins that I used to partition mutations.

For now, I’ll combine these into a single cluster. I’ll explore this further with bridging reads later in the analysis.

updated.haplotype.df = updated.haplotype.df %>%
  mutate(Haplotype = case_when(Haplotype %in% c("cluster 14") ~ "cluster 10",
                               TRUE ~ Haplotype))

I’ll write out these newly phased clusters of mutations for further analysis.

## [1] "Writing the results to results/variants/clustered_variants.csv"

Visualize All Mutations

At this stage in the analysis, we’ve phased mutations into rough haplotypes. We’ve also assigned some mutations to the major ‘genotypes’ (Genome 1 and Genome 2) that were missed in the previous notebook due to low coverage.

Let’s visualize all of these mutations in each tissue specimen. Mutations are annotated by whether they’re Genome 1, Genome 2, or a ‘Driver’ mutations. These ‘Drivers’ are mutations with impacts similar to mutations that have been observed in previous cases of SSPE.

What are the SNPs we’re calling ‘Genome 1 (G1)’. These will be different than the ‘Candidate Genome 1’ SNPs in Figure 2, because these were identified by their correlation in frequency between all tissue specimens.

What about Genome-1-1?

In the previous notebook, we noticed a set of mutations that made up a haplotype that arises on the background of Genome 1, but it mostly missing from the Frontal Cortex 2 sample. We called this haplotype Genome-1-1, and it’s very significant because it’s indicative some early divergence in the frontal cortex. Let’s take a look at this haplotype in every tissue.